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This  work  discusses  a  Bayesian  approach  to  approximating  the  distribution  of 
parameters  governing  nonlinear  structural  systems.  Specifically,  we  use  a  Markov 
Chain  Monte  Carlo  method  for  sampling  the  posterior  parameter  distributions  thus 
producing  both  point  and  interval  estimates  for  parameters.  The  method  is  first  used  to 
identify  both  linear  and  nonlinear  parameters  in  a  multiple  degree-of-freedom 
structural  systems  using  free-decay  vibrations.  The  approach  is  then  applied  to  the 
problem  of  identifying  the  location,  size,  and  depth  of  delamination  in  a  model 
composite  beam.  The  influence  of  additive  Gaussian  noise  on  the  response  data  is 
explored  with  respect  to  the  quality  of  the  resulting  parameter  estimates. 
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1.  Introduction 

System  identification  can  be  loosely  defined  as  the  process  of  estimating  parameters  associated  with  a  specified  model 
(or  models)  from  acquired  data.  There  are  two  main  schools  of  thought  in  estimation  problems:  the  frequentist  approach, 
often  based  on  the  method  of  maximum  likelihood  (ML),  and  the  Bayesian  approach.  Both  methods  seek  to  provide  the 
best  possible  parameter  estimates  in  the  face  of  the  inevitable  uncertainty  (e.g.  measurement  error)  present  in  the 
observed  data.  In  fact  the  likelihood  function,  describing  the  joint  distribution  of  the  data  given  the  model  and  model 
parameters,  is  the  central  ingredient  in  both  approaches  to  estimation.  However,  the  two  approaches  are  fundamentally 
different  in  how  they  treat  model  parameters  resulting  in  different  approaches  to  inference. 

Historically,  researchers  working  on  nonlinear  system  identification  problems  have  tended  to  focus  on  ML  methods.  In  a 
recent  review  paper  Kerschen  et  al.  [1,  see  Sections  6  and  7  and  the  many  references  contained  therein],  cover  numerous 
available  techniques  for  nonlinear  parameter  estimation.  These  approaches  are  based  on  time  domain  (e.g.  [2,3]), 
frequency  domain  (e.g.  [4-6])  or  higher-order  spectral  analysis  [7,8].  Each  of  these  techniques  takes  the  best  estimate  to  be 
the  one  that  minimizes  the  mean  square  error  between  data  and  model.  Although  not  explicitly  stated,  this  choice  of  cost 
function  yields  ML  estimates  for  the  parameters,  provided  that  the  uncertainty  in  the  model  is  taken  as  additive,  iid,  jointly 
Gaussian  noise.  The  specific  pros  and  cons  associated  with  these  methods  are  highlighted  in  [1]  and  are  therefore  not 
discussed  here.  Rather,  this  work  is  focused  on  the  conceptually  different  Bayesian  alternative  to  nonlinear  parameter 
estimation  problems. 
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Perhaps  the  chief  benefit  of  the  Bayesian  worldview  is  that  it  treats  the  parameters  we  seek  to  estimate,  denoted  by  the 
vector  0,  as  random  variables  with  joint  distribution  p(0).  This  allows  us  to  make  probabilistic  statements  about  our 
estimated  parameters.  For  example,  let  us  say  our  goal  is  to  collect  data  from  a  structure  and  estimate  the  crack  length, 
denoted  “  0i”  .  The  Bayesian  framework  allows  us  to  make  statements  such  as  the  following: 


There  is  a  95  percent  chance  that  the  crack  length  lies  in  the  interval  0.1  <  <  0.5  [cm]. 


There  is  no  way  to  make  an  analogous  statement  using  the  frequentist  view  of  the  world  because  the  crack  length 
would  not  be  treated  as  a  random  variable.  Instead  we  would  have  to  make  statements  about  the  repeatability  of  our 
estimation  procedure  i.e.  the  statistical  machinery  that  produced  the  interval  estimate.  Thus  we  might  be  able  to  make  the 
statement  “The  methods  we  have  used  produce  random  intervals  which  may  or  may  not  include  the  unknown  parameter; 
they  succeed  in  95  percent  of  similar  circumstances,  but  we  can’t  say  anything  with  certainty  about  [0.1, 0.5]’’  Although  in 
some  cases  the  two  methods  will  produce  the  same  estimate  and  interval  the  interpretation  of  the  interval  is  different,  the 
Bayesian  interval  being  a  probability  statement  about  the  parameter.  Thus,  our  main  reason  for  gravitating  toward  the 
Bayesian  viewpoint  is  that  it  provides  a  direct  estimate  of  the  information  we  are  after  as  opposed  to  related  (but  different) 
information.  Additional  reasons  for  the  Bayesian  approach  are  highlighted  in  Chapter  1  of  [9].  First,  the  approach  is  more 
amenable  to  more  complex  data  and/or  models  than  is  the  frequentist  (ML)  approach  (including  cases  where  a  frequentist 
approach  does  not  exist).  Secondly,  the  approach  described  here  does  not  require  asymptotics  in  approximating  the 
confidence  intervals  (as  does  the  frequentist  approach  of  using  the  Fisher  Information  Matrix  to  bound  the  variance  of  a 
parameter  estimate,  see  e.g.,  [10]). 

This  is  certainly  not  the  first  paper  to  explore  Bayesian  methods  in  structural  dynamics.  A  general  Bayesian  approach  to 
structural  dynamics  problems  was  put  forth  in  the  ’90s  with  works  by  Beck  and  Katafygiotis  [11,12].  These,  and  subsequent 
works  (see  e.g.,  [13]),  relied  on  asymptotics  in  order  to  solve  for  the  marginal  parameter  distributions  which,  as  will  be 
shown,  often  take  the  form  of  high-dimensional  integrals.  At  around  the  same  time  Sohn  and  Law  [14]  also  developed  a 
Bayesian  approach  to  the  structural  identification  problem.  Their  approach,  also  described  in  subsequent  papers  [15,16], 
circumvented  the  issue  of  solving  for  the  desired  marginal  parameter  distributions  by  using  a  so-called  “branch-and- 
bound’’  strategy  to  instead  solve  for  the  most  likely  damage  hypothesis  in  a  damage  detection  application. 

Rather  than  abandon  the  goal  of  parameter  estimation  or  rely  on  asymptotics,  one  can  make  use  of  a  powerful  approach 
for  sampling  the  marginal  distributions  without  having  to  perform  the  integration.  Such  an  approach  was  first  proposed  by 
Metropolis  et  al.  [1 7]  in  the  early  1950’s  for  solving  high-dimensional  integrals  in  particle  physics.  A  later  work  by  Hastings 
[18]  extended  this  general  approach,  resulted  in  what  is  now  known  as  the  Markov  Chain  Monte  Carlo  (MCMC)  approach  to 
approximating  probability  distributions.  MCMC  methods  have  since  become  extremely  popular  in  implementing  Bayesian 
estimation  in  a  number  of  fields  ranging  from  ecology  [9]  to  genetics  [19]  and  have  recently  seen  increased  use  in 
engineering  applications.  Beck  et  al.  appear  to  have  pioneered  the  use  of  MCMC  methods  in  structural  dynamics  in  Beck 
and  Au  [20]  and  more  recently  in  Cheung  and  Beck  [21].  Additional  work  by  Glaser  et  al.  [22]  illustrated  the  approach  in 
detecting  stiffness  reduction  in  beams  using  static  measurements.  This  method  has  also  been  used  to  estimate  failure 
probabilities  in  structural  reliability  problems  as  part  of  the  “Subset  Simulation’’  approach  of  Au  et  al.  [23].  In  related 
works,  Zhang  and  Cho  [24]  used  the  MCMC  approach  to  help  design  an  evolutionary  algorithm  for  performing  system 
identification,  while  Kerschen  et  al.  [25]  used  the  MCMC  approach  to  select  among  competing  models  for  describing  the 
dynamics  of  a  nonlinear  mechanical  system. 

Our  goal  in  this  work  is  to  focus  on  the  use  of  the  MCMC  approach  to  Bayesian  parameter  estimation  in  nonlinear 
systems  using  free-vibration,  time-series  data.  The  approach  is  certainly  more  general  and  could  be  applied  to  forced 
structures  as  well.  However,  in  practice  the  forcing  function  is  not  always  obtainable  while  for  the  free-decay  problem  we 
simply  have  to  include  the  initial  conditions  as  random  variables  to  be  predicted.  A  different  approach  that  does  not 
require  input  data  would  be  to  perform  the  analysis  using  estimated  frequencies  as  a  means  of  comparing  model  to  data  as 
was  done  in  [13].  However  here  the  practitioner  has  the  additional  step  of  estimating  the  frequencies  from  observed  data. 
This  is  compounded  by  the  task  of  determining  the  analytical  model  frequencies  which  for  nonlinear  systems  can  be 
extremely  challenging.  Thus,  free-decay  response  data  provide  for  a  direct  model-to-data  comparison  and  are  easily 
obtainable  in  experiment.  Another  practical  advantage  of  this  general  approach  to  system  identification  is  that  one  does 
not  need  to  measure  time-series  data  from  each  of  the  structural  degrees-of-freedom  (DOF)  in  order  to  estimate  the 
associated  parameters.  For  example,  a  single  time-series  response  from  one  of  the  DOFs  can,  in  some  cases,  be  used  to 
estimate  model  parameters  associated  with  DOFs  not  observed. 

In  this  work  the  approach  is  first  used  to  estimate  the  parameters  associated  with  a  two  degree-of-freedom  nonlinear 
structural  system.  The  relationship  between  the  quality  of  the  resulting  estimates  and  the  signal-to-noise  ratio  is  explored. 
We  then  turn  our  attention  to  the  difficult  problem  of  estimating  and  tracking  delamination  growth  in  a  composite  beam 
model.  This  model  was  recently  developed  in  [26]  where  it  was  shown  to  accurately  capture  the  localized  buckling  that 
occurs  due  to  the  separation  of  the  laminates.  Subsequent  work  by  the  authors  [27]  focused  on  detection  of  the 
delamination  using  a  higher-order  spectral  analysis.  Here,  the  focus  is  on  identifying  the  damage  parameters  using  only 
free  vibration  data,  a  task  to  which  the  Bayesian  approach  using  MCMC  to  sample  the  desired  parameter  distributions  is 
well-suited. 
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2.  Bayesian  approach  to  structural  parameter  estimation 


In  structural  dynamics  the  uncertainty  that  gives  rise  to  the  likelihood  function  is  typically  the  result  of  measurement 
error  on  the  observed  data.  Assume  that  we  have  measured  a  structural  response  at  a  specific  location  giving  the  N-point 
time-series 


Zn=yn  +  ^„,  n=l...N  (1) 

where  yn  is  the  noise-free  response  and  is  a  sequence  of  independent,  identically  distributed  (iid)  samples  drawn  from  a 
zero-mean  Gaussian  distribution  with  variance  In  this  case  we  can  write  the  likelihood 


Pl(z|0,  al)  = 


1 

(2nal)^^^ 


exp 


1 

2al 


N 


(z„-y„(0))2 

n  =  1 


(2) 


where  we  have  specified  the  dependence  on  the  model  parameters  through  the  model  response  yn(0).  The  likelihood  function 
Pl  (z|0,  (Tq)  is  a  probabilistic  statement  (model)  about  the  distribution  of  the  observed  data  z  =  Zn  n  =  l  ...N  given  a  model 
response,  determined  by  the  parameter  vector  0.  The  method  of  maximum  likelihood  estimates  the  model  parameters  to  be 
those  that  maximize  Pl(z|0,  cr^).  Some  important  and  useful  properties  of  MLEs  is  that,  under  regularity  conditions,  they 
are  asymptotically  unbiased  and  possess  the  minimum  possible  variance  in  the  resulting  parameter  estimates.  However,  the 
approach  suffers  from  the  practical  and  conceptual  drawbacks  already  mentioned  in  the  introduction. 

In  the  event  that  data  from  M  structural  locations  are  recorded  the  likelihood  function  can  be  easily  altered  as 


Pl(z|0,  al)  = 


1 


1 

2al 


M  N 

Y  Y(^r-yT\m^ 


m  =  1  n  =  1 


(3) 


where  the  superscript  (m)  will  henceforth  be  used  to  denote  the  degree-of-freedom  being  observed.  In  this  case  we  have 
made  two  important  but  realistic  assumptions.  First,  we  have  assumed  that  the  sensor  noise  is  uncorrelated.  Secondly  we 
have  assumed  that  the  sensor  noise  has  the  same  variance  across  all  sensors.  Neither  assumption  is  critical  to  the  proposed 
approach — for  example  if  the  sensor  noise  was  assumed  different  for  each  sensor  we  would  simply  identify  each  noise 
variance  separately  in  the  algorithm  to  be  described  in  Section  3.  This  same  likelihood  function  was  used  in  [21].  It  should 
be  mentioned  that  the  approach  does  not  require  all  degrees-of-freedom  to  be  measured,  that  is  to  say  M  <  DOF.  Of  course 
the  quality  of  the  estimates  typically  improves  with  data  from  multiple  measurement  points,  but  is  not  strictly  necessary. 
Finally,  we  should  also  point  out  here  that  this  framework  can  accommodate  multiple  models  Mi  i  =  1 . . .  Nm  such  that 
the  likelihood  function  becomes  pL(z\d,aQ,Mi)  (i.e.  the  likelihood  is  conditional  on  the  specified  model).  However,  in  this 
work  we  will  only  consider  a  single  model  and  drop  the  model  notation  Mi.  Additionally,  in  the  following  general 
discussion  we  consider  the  parameter  vector  to  include  the  noise  variance  for  notational  convenience. 

Bayesian  analysis  relies  on  Bayes  Theorem,  relating  the  prior  parameter  probability  density  function  (PDF)  Pn(^)  to  the 
posterior  p(0|z)  via 

p(0|z)  =  Pi(z|0)P;,(0)/PD(z).  (4) 


The  key  ingredient  relating  the  prior  to  posterior  distributions  is  the  likelihood  function.  Eq.  (4)  provides  us  a  simple  means 
of  relating  prior  information  to  the  parameter  distributions  we  want.  The  term  in  the  denominator,  Pd(z),  is  a  normalizing 
constant  that  can  be  ignored  in  the  following  development.  The  joint  parameter  distribution  contains  the  desired  marginal 
distributions  of  each  of  the  P  parameters 

P(0j\z)=  f  ^  p(0|z)d0_joc  f  pL(z|0)P;t(0)d0_j  (5) 

where  the  notation  d0_j  denotes  the  multidimensional  integral  over  all  parameters  other  than  Oj.  Thus,  we  require  a 
means  of  performing  a  potentially  high-dimensional  integral  involving  a  likelihood  function  for  which  we  will  often  have  a 
very  complicated  expression.  For  example,  an  analytical  expression  for  yn  for  a  nonlinear  system  response  is  a  very 
challenging  problem  for  even  a  single  DOF  system  possessing  simple  nonlinearities.  Fortunately  there  exists  a  convenient 
numerical  approach  to  sampling  from  the  marginal  parameter  distributions. 


3.  Markov-chain  Monte-Carlo  methods 

The  term  “Monte  Carlo  methods”  is  used  to  describe  simulation  techniques  for  investigating  probability  distributions. 
These  techniques  can  be  highly  efficient,  especially  when  independent  samples  can  be  generated.  Unfortunately,  posterior 
distributions  used  in  Bayesian  inference  are  often  complicated,  making  it  difficult  to  draw  independent  samples. 
Nevertheless  it  is  often  easy  to  draw  a  dependent  sequence  of  samples  representing  posterior  distributions.  Over  the  last 
60  years,  stochastic  algorithms  have  been  developed  which  sample  new  values  using  rules  determined  by  a  fixed  number 
of  previous  observations.  The  result  is  a  Markov  chain;  this  form  of  simulation  is  known  as  Markov  chain  Monte  Carlo 
(MCMC).  The  magic  of  MCMC  is  in  producing  algorithms  for  which  the  resultant  Markov  chains  have  stationary 
distributions  equal  to  the  distribution  we  wish  to  sample. 
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3.1.  Gibbs  sampling 

Gibbs  sampling  is  a  form  of  MCMC  used  to  sample  multivariate  distributions.  Consider  a  vector  valued  parameter 
d  =  (0-i,02,. . .  ,9py,  data  z,  prior  P7r(0),  and  likelihood  Pl(z|0).  Define  0_j  as  the  vector  obtained  by  deleting  the  jth 
component  of  0,  viz., 

0_j  =  (01 , 02,  •  •  •  ,  0j-l ,  0j+ 1 , .  .  .  ,  Op)'. 

Gibbs  sampling  of  the  posterior  distribution 

P(0|z)ocPl(z|0)P;,(0) 

produces  a  sequence  (Markov  Chain)  {0®}i  =  1 ...  ChainLength  with  components  0-*^  generated  sequentially,  for  j  =  l,2,...,P. 
The  value  0j'^  is  sampled  from  its  full  conditional  distribution  p(0j'^|z,0_j  =  which  for  simplicity  is  denoted  as  p(0®|  •). 
The  full  conditional  distribution  of  Oj  can  be  thought  of  as  the  posterior  distribution  of  0j,  if  it  happened  that  the  other  P  -1 
parameters  were  known,  and  equal  to  their  present  values,  i.e.,  if  6^  =  0^^*“^^  for  k^j.  Like  the  posterior  distribution  p(0|z), 
the  full  conditional  distribution  is  proportional  to  the  product  of  likelihood  and  prior, 

p(0f\  •)ocpi(z|0)p^(0): 

the  difference  is  that  p(0®|  •)  is  a  function  of  Oj  alone,  rather  than  of  all  P  variables  comprising  0. 

To  implement  Gibbs  sampling,  then,  we  only  need  to  be  able  to  sample  full  conditional  distributions.  The  ease  or 
difficulty  of  this  endeavor  depends  on  the  form  of  the  likelihood  function  and  the  prior.  Among  the  easy  cases  is  the 
circumstance  of  conditional  conjugacy,  when  full  conditional  and  prior  are  both  identifiable  as  members  of  the  same 
parametric  family. 

3.2.  Conditional  conjugacy 

To  illustrate  the  notion  of  conjugacy,  consider  z,  the  number  of  successes  in  N  independent  Bernoulli  trials,  with  success 
parameter  0.  A  likelihood  function  for  0  is 

Pl(Z|0)  =  0^(l-0)'^-". 

A  beta  prior  for  0  on  [0,1]  has  density  function 

p^(0)oc0“-’(l-0/-\ 

for  some  values  a,  >  0.  The  posterior  distribution  is  proportional  to  the  product  of  likelihood  and  prior,  hence  is  clearly  a 
member  of  the  beta  family  of  distributions,  with  parameters  a'  =  (x-\-z  and  jf'  =  jf+N-z.  Beta  priors  are  conjugate  to  the 
binomial  likelihood. 

Suppose  that  the  parameters  and  0  are  unknown  in  model  (2)  or  (3).  For  reasons  that  will  become  clear 
subsequently,  we  reparameterize  the  model  using  the  precision  parameter  t  =  1  /cr^,  and  show  how  conjugacy  can  be  used 
in  describing  a  full  conditional  distribution  for  t.  Given  independent  priors  on  Pni^:)  and  Pn(^),  the  posterior  distribution  is 

p(0,  T|  z) cc  y  Q(z,  0))p„(T)P;t(0), 

where 


M  N 

a(z,  0)  =  ^  ^  (6) 

i  =  1  n  =  1 

is  the  sum  of  squares  in  the  likelihood.  It  follows  that  the  full  conditional  distribution  for  t  is 

p(T|  z,  0)  oc  y  Q(z,  0)) p,(T). 

The  gamma  family  of  distributions  is  described  by  density  functions  p(x)ocx°^“^exp(-j5x)  for  x>0,  indexed  by 
parameters  a,p>0.  Thus  if  t  has  a  gamma  prior  with  parameters  a  and  jf,  the  full  conditional  distribution  for  t  is  also  in 
the  gamma  family,  with  parameters  oc'  =  a+MN/2  and  P'  =  j5  +  Q.(z,0)/2.  The  gamma  prior  is  conditionally  conjugate  for  t. 
We  may  therefore  directly  sample  t  (hence  cr^)  at  each  step  in  the  Markov  chain  without  having  to  resort  to  the  more 
computationally  intensive  sampling  approach  described  in  the  next  section.  In  this  work  we  use  the  so-called  “diffuse” 
(uninformative)  prior,  obtained  by  setting  a  =  \,p  =  0.  Conditional  conjugacy  simplifies  Gibbs  sampling  by  allowing 
reference  to  standard  probability  distributions  (like  the  beta  or  gamma)  readily  available  in  standard  software  packages. 

3.3.  Metropolis-Hastings  algorithm 

Unfortunately,  conjugate  or  conditionally  conjugate  priors  are  not  always  available.  In  such  cases,  techniques  such  as 
rejection  sampling  can  be  used  for  sampling  the  full  conditional  distribution.  One  may  also  resort  to  the  omnibus 
Metropolis-Hastings  algorithm,  a  simple  implementation  of  which  we  now  describe. 
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Suppose  that  we  wish  to  draw  samples  from  a  specific  distribution /(0)  (in  our  case  f(0)  might  be  the  left-hand-side  of 
Eq.  5,  p(9j)).  The  Metropolis-Hastings  algorithm  generates  a  sequence  {d®}  by  a  2-step  procedure.  At  stage  i,  a  candidate 
value  0*  is  sampled  based  on  the  current  value  it  is  sampled  from  a  distribution  function  A  Bernoulli  trial  is 

performed,  with  success  probability  f  =  min{r,  1},  where 

if  the  result  of  the  trial  is  success,  we  set  6^'^  =  9*',  otherwise,  we  set  9^'^  =  9^'~^\  The  alternatives  are  often  described  as 
“moving”  or  “staying  at  the  current  value”,  but  it  is  important  to  note  that  in  the  latter  case,  the  same  value  is  recorded 
twice  (once  as  once  as 

It  is  often  easy  to  describe  a  candidate  generating  distribution  which  is  symmetric  in  its  arguments,  hence  cancels  out  of 
the  numerator  and  denominator  of  r.  For  example,  if  0*  is  uniformly  distributed  on  an  interval  of  length  2A  centered  on 
then  its  density  function  is 


(7) 


so  that  g(0<''"’’|r)  =g(0<''"’’|r). 

Suppose  that  we  use  this  uniform  candidate  generator,  with  A=l,  and  apply  the  technique  to  generating  samples  from 
the  standard  normal  distribution /(d)  =  exp(-d^/2)/V^.  Then  r  simplifies  to 


r  =  exp 


Two  features  of  this  remarkable  algorithm  deserve  comment.  First,  that  the  first-order  Markov  chain  it  produces 
will  have  stationary  distribution  equal  to  f(0),  provided  only  that  the  candidate  generating  distribution  and  the 
starting  value  0^^^  allow  the  chain  to  reach  all  values  t  for  which  f(6)  >  0.  Thus  in  generating  standard  normal  variates 
with  the  uniform  candidate  generator,  it  does  not  matter  if  we  start  the  chain  at  =  25,  even  though  such 
values  are  exceedingly  rare  for  a  standard  normal  variate.  Most  of  the  first  50  or  so  moves  toward  zero  will  all  be 
accepted,  while  most  of  the  first  50  or  so  moves  away  from  zero  will  all  be  rejected,  and  the  chain  will  rush  in  towards 
the  range  (-3,3)  typical  of  standard  normal  variables.  These  initial  values  are  discarded  in  MH  sampling,  described  as  a 
“burn-in”. 

The  second  important  feature  of  this  algorithm  is  that  its  implementation  depends  on  the  desired  distribution  only 
through  the  ratio  in  r.  Thus  in  calculating  r  for  the  standard  normal  distribution,  the  constant  terms  V2n  cancel.  This  is  an 
enormous  boon  in  calculating  full  conditional  distributions,  which  are  defined  implicitly  as  the  distribution  proportional  to 
product  of  likelihood  and  prior;  the  normalizing  constant  in  Eq.  (4),  which  makes  the  distribution  integrate  to  1,  need  not 
be  calculated. 

Practical  implementation  of  the  MCMC  technique  can  require  some  tuning.  Large  values  for  the  “tuning  parameter”  A  (Eq.  (7)) 
decreases  the  likelihood  of  accepting  a  new  value,  thus  the  samples  in  the  chain  tend  to  be  highly  correlated.  Too  small  a  step, 
however,  and  the  algorithm  can  take  a  prohibitively  long  time  to  converge.  The  tradeoff  between  excessive  computation  time 
(small  steps)  and  not  generating  independent  samples  (big  steps)  is  well-known.  Here  we  use  a  simple  approach  whereby  the 
tuning  parameter  is  adjusted  on-the-fly  in  order  to  achieve  an  appropriate  acceptance  rate  of  between  30  percent  and  50  percent. 
This  range  of  acceptance  probabilities  has  been  demonstrated  to  produce  Markov  chains  with  low  auto-correlation  (good  mixing) 
[9].  During  the  burn-in  period  we  simply  divide  A  by  the  constant  value  1.007  after  each  rejection  and  multiply  by  1.01  after  each 
acceptance.  Thus,  an  acceptance  causes  us  to  expand  our  parameter  search  while  rejection  results  in  smaller  “kicks”  to  the 
previous  value  in  the  chain.  The  asymmetry  in  the  constants  causes  a  slight  bias  in  the  favor  of  rejection.  This  simple  approach  is 
quite  effective  at  producing  acceptance  rates  of  40  percent,  a  good  target  for  the  MCMC  algorithm.  Details  of  this  approach  are 
discussed  in  [9]. 


Algorithm  1.  The  MCMC  algorithm  using  Metropolis-Hastings  with  Gibbs  sampling  for  Gaussian  likelihood  and  Uniform 
“transition”  distribution  g(0*,  dj(i-l))  =  1  /(2Aj) 

Task:  Generate  the  posterior  parameter  distributions  p(9j)  for  j  =  1 . .  .P  given  the  model,  the  observed  data  z  =  Zn  n^l  ...N,  and  prior  parameter 
distributions  Pnjidj)-  Also  estimate  the  noise  variance  Op+i  = 

Initialization:  Initialize  chain  i=0 

Initial  guesses  for  parameter  values  (chosen  from  the  priors)  0(0)  =  0j{O)  j  =  1 , .  .P 
Initial  values  for  tuning  parameters  Aj  ...P 
Sample  initial  variance  from  inverse  Gamma  distribution 
Op^  1  (0)  =  lG(N/2  + 1 , 0.5Q(z,  0(0)) 

Set  number  of  burn-in  iterations  B 
Main  iteration: increment  i  by  1  and  apply 
For  each  parameter  j 

Generate  candidate  6*  =  6j(i-l)  +  2Aj  x  U(-l,  1)  where  U{a,b)  is  a  uniformly  distributed  number  on  [a,b] 

Compute  r  =  exp((-0.5/ep+i(i-l))  x  ((2(z,0*)-(2(z,0(i-l))))  *ppj(0*)/ppj(0(i-l))  where  0*  =(0i(O, M'-l))  and 
0(!  - 1 )  ^  (S,  (i),  •  •  ■ ,  Sj- ,  (0,  e/i- 1 ),  + 1  (i- 1 ),  ■  ■  ■ ,  ep(!- 1 )) 
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If  U(0, 1)  <  r  keep  the  new  value,  i.e.  set  9j(i)  =  6*,  and  adjust  tuning  parameter  Aj  ^Aj  x  1.01 
Else  reject  the  new  value,  keeping  9j(i)  =  9j(i-\)  and  adjust  tuning  parameter  Aj=Ajl\. 007. 
Directly  sample  the  variance  posterior  9p+ 1  (I)  =  IG(N/2  + 1 , 0.5Q(z,  0(i)) 

After  i  >  B  iterations,  cease  adjusting  the  Aj  and  record  subsequent 
values  9j(i)  as  members  of  p(dj). 


The  algorithm  for  executing  this  specific  implementation  of  MCMC  is  given  in  Algorithm  1.  This  pseudo-code  implicitly 
assumes  (1)  a  diffuse  Gamma  prior  on  the  precision  parameter  l/cr^,  (2)  a  Gaussian  likelihood  function,  governed  by  the 
sum-squared  error  between  data  and  model  and  (3)  a  uniform  transition  distribution  (Eq.  (7)).  The  code  is  quite  simple  to 
implement,  the  costly  step  being  the  determination  of  the  sum-squared  error  function  Q.(-,  )  (Eq.  (6)).  Eor  nonlinear 
systems,  for  which  no  analytical  solution  exists,  generating  the  model  output  yn  for  a  given  parameter  vector  requires 
numerically  integrating  the  equations  of  motion.  Also  note  that  in  the  above  described  implementation,  as  the  Gibbs 
sampler  moves  through  the  parameter  vector  the  newly  updated  values  are  used.  That  is  to  say,  for  j=2  we  use  6\{i)  rather 
than  0i(i-l)  as  the  first  element  of  0(i-l).  Einally,  many  statistical  packages  provide  routines  that  generate  samples  from  a 
Gamma  distribution  but  not  from  an  inverse  Gamma  distribution.  Eortunately  there  is  a  simple  relationship  between  the 
Gamma  G(-,  •)  and  inverse  Gamma  IG(-,  •)  distributions 

IG(a,  b)  =  1 /G(a,  1 /b) 

so  that  for  our  problem  we  may  sample  the  variance  by  drawing 

1  (0  ~  1 .0/[G(N/2  -E 1 , 2/(2(z,  0(0)]. 

There  are  obviously  many  variations  on  the  above  described  algorithm,  however  Algorithm  1  presents  a  simple  but  useful 
implementation  that  is  straightforward  to  code  in  software. 


4.  Example  1:  2-DOF  nonlinear  spring-mass-damper  system 


In  order  to  illustrate  the  above  described  identification  procedure,  consider  the  two  degree-of-freedom  (DOE)  system 
described  by  the  second  order,  nonlinear,  ordinary  differential  equations 


[M]yt + [Qyt + [/qyt +g(yt,  y,)  =  0 


(8) 


where 


[M]  = 


mi 

0 


0 

m2 


[C]  = 


C1+C2  -C2 

-C2  C2 


[K]  = 


/<1-E/<2  -/<2 

-/<2  /<2 


are  constant  coefficient  mass,  damping  and  stiffness  matrices  respectively.  The  nonlinear  function  g(  )  provides  quadratic 
coupling  between  masses.  Here  we  consider  a  quadratic  restoring  force  between  masses  1  and  2  so  that 


g(yt,yf)  = 


-k„o.(y?>-y?V 

k„on(y?^-y?b^ 


where  knon  is  the  nonlinear  stiffness  coefficient.  We  consider  as  the  observed  data  N=512  points  of  a  noise-contaminated, 
free-decay  response,  obtained  by  numerically  integrating  Eq.  (8)  with  a  time-step  of  zlt  =  0.01s  so  that  discrete  and 
continuous  time  are  related  via  t  =  nAt.  The  initial  conditions  were  set  to  the  values  y[/^  =  0.0m,  =  O.Om/s, 

yQ^^  =  0.15  m,  and  =yQ^^  =  O.Om/s.  The  parameter  values  used  in  generating  the  response  data  were  /<:i  =  2000N/m, 
/<2=1500N/m,  Cl  =C2  =  6Ns/m,  and  /<non=l 0,000 N/m^.  The  mass  parameters  were  fixed  to  the  values  mi  =  m2  =  l  kg  and 
were  assumed  known  (measured)  at  the  outset.  The  variance  of  the  additive,  Gaussian  distributed  noise  is  defined  by  the 
signal-to-noise  ratio  SNR  =  aj /a\  so  that  the  observed  data  are  given  by 

where  each  is  an  iid  Gaussian  random  variable  and  aj  is  the  variance  of  the  true  underlying  signal.  We  are  assuming  in 
this  example  that  only  a  single  time-series  measurement  is  available  in  forming  the  likelihood  (i.e.  M=1  in  Eq.  (3)).  Eor  this 
example  the  additive  noise  level  was  set  at  -  lOdB  down,  i.e.  the  signal  variance  was  10  x  greater  than  the  noise  variance. 
The  job  of  the  algorithm  is  to  identify  Ci,C2,/<i,/<2  as  well  as  the  nonlinear  parameter  know  the  initial  conditions  y^'^\ 
and  the  noise  variance  a^. 

Eig.  1  shows  the  resulting  stationary  parameter  distributions  for  each  of  the  identified  quantities.  Each  of  the  model 
parameters  were  identified  using  the  above-described  Metropolis-Hastings  algorithm  with  Gibbs  sampling  and  assuming 
Uniform  priors.  The  exception  was  the  variance  a^  which  could  be  sample  directly  by  assuming  a  conjugate  Gamma  prior 
as  described  in  Section  3.2.  In  other  words,  at  each  step  in  the  chain  we  choose  ~  l/G(MN/2,-^^)  where  G(  )  is  the 
Gamma  distribution.  In  building  the  Markov  Chains  we  used  a  burn  in  of  1 50,000  iterations  and  retained  the  subsequent 
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Fig.  1.  True  model  output  and  corrupted  (observed)  data  (a)  followed  by  the  posterior  distributions  of,  respectively  (b-k),  /<i,  /C2,  Ci,  C2,  y^\  y^\  v^\ 

crj,  and  knon- 
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Fig.  2.  Progression  of  the  95  percent  confidence  interval  associated  with  the  parameter  estimates  as  a  function  of  signal-to-noise  ratio:  (a)  linear  stiffness, 
/<!  (b)  linear  stiffness,  ^2  (c)  linear  damping,  Ci  (d)  linear  damping,  C2  (e)  nonlinear  stiffness,  knon  and  (f)  noise  variance,  As  expected,  a  higher  SNR  tends 
to  result  in  better  estimates. 
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50,000  iterations  as  the  stationary  posterior  distributions  of  interest.  Fig.  1  shows  the  noise-corrupted  time-series  data 
used  in  the  identification  along  with  all  identified  parameters.  The  algorithm  is  clearly  able  to  identify  both  the  system 
parameters  and  the  initial  conditions  using  only  a  short,  noisy  free-decay  response  as  the  observed  data.  Certain 
parameters  are  identified  with  more  confidence  than  others.  Notoriously  difficult  to  obtain  damping  estimates  have  the 
largest  percentage  variance  of  all  parameters  while  the  stiffness  estimates  are  better  defined.  The  nonlinearity  parameter  is 
similarly  estimated  to  lie  within  a  fairly  narrow  confidence  band. 

As  expected,  the  distributions  tend  to  narrow  as  the  signal-to-noise  ratio  is  increased.  Fig.  2  shows  the  progression  of 
the  95  percent  confidence  intervals,  defined  as  the  central  95  percent  of  the  50,000  values  that  comprise  the  stationary 
Markov  chain,  as  a  function  of  SNR.  Using  this  approach  each  of  the  parameters  is  correctly  identified  with  their  associated 
confidence  intervals.  There  is  no  need  to  base  the  confidence  intervals  on  an  assumed  Gaussian  distribution  as  is  often 
done  in  parameter  estimation  (e.g.  assume  +  3  standard  deviations  from  the  mean  as  95  percent  confidence).  The  MCMC 
algorithm  provides  an  approximation  of  the  entire  posterior  PDF,  regardless  of  its  underlying  form.  The  above  results  have 
demonstrated  that  we  can  successfully  identify  both  linear  and  nonlinear  system  parameters  in  a  multi  DOF  structure 
using  only  measurements  of  noisy,  free-decay  response  data.  The  next  example  involves  the  more  difficult  problem  of 
detecting  delamination  in  composite  structures. 


5.  Delamination  identification  in  a  composite  beam 


As  a  second  example,  we  seek  to  identify  the  location,  extent,  and  depth  of  delamination  in  a  composite  beam  structure. 
The  dynamic  beam  model  was  derived  previously  in  [27]  and  is  shown  schematically  in  Fig.  3.  This  model  is  low 
dimensional  (only  three  independent  coordinates  need  to  be  specified),  yet  was  shown  experimentally  to  accurately 
capture  the  localized  buckling  that  occurs  due  to  the  presence  of  the  delamination  under  static  loading  [26].  The  global 
beam  motion  is  assumed  to  be  dictated  by  the  first  mode  of  the  response,  thus  the  global  displacement  assumes  the  form 

where 


•AiW 


©■-5© 


31 


is  a  normalized  shape  function  describing  the  vertical  beam  deflection  at  any  point  x,  measured  from  the  left  end  of  the 
beam.  The  other  two  coordinates,  describing  the  time-dependent  motion  of  the  upper  (region  2)  and  lower  (region  3) 
laminates  respectively,  are  assumed  to  be  of  the  form 

y2(X2,t)=y^(X2+Xa,t)+q2(t)'P2(X2)+^(1-a)h 


y3(X3,t)  =yi(X3  +Xa,  t)  +  q3(t)'P3(X3)-  y 


<^2,3  =  1-COS^ 


Xb-Xa 


(9) 


where  the  shape  functions  ^2,3  describe  the  deflected  shape  of  regions  2  and  3  in  Fig.  3.  The  constant  terms  in  Eq.  (9) 
denote  the  neutral  axis  offsets,  measured  from  the  global  neutral  axis,  of  each  of  the  laminates.  Substituting  these 
expressions  into  Lagrange’s  equation  yields  a  set  of  three  coupled,  nonlinear  differential  equations  in  terms  of  the 
time-dependent  vector  q(t)  =  (qi(t),  q2(t),  Qsit)) 

[M]qt  +  [C]qr  +  [Ki]qr  +  [K^lqi  (Oqt  +  +  [Kc]q^  =0.  (10) 


Fig.  3.  Schematic  of  the  dynamic  delaminated  beam  model.  Region  1  is  simply  modeled  as  a  linear  cantilevered  beam  whose  motion  is  governed  by  the 
first  mode  of  vibration.  Regions  2  and  3  are  modeled  as  nonlinear  beams  where  axial  stretching  is  permitted. 
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The  equations  of  motion  are  therefore  described  in  terms  of  the  material  properties  of  the  beam,  the  dimensions  of  the 
beam,  and  the  parameters  associated  with  the  damage.  These  damage  parameters  are  Xa,  the  delamination  starting 
location,  X/,,  the  delamination  end  point,  and  a,  the  delamination  depth  specified  as  a  fraction  of  the  overall  beam  thickness 
h.  Additionally,  we  have  chosen  to  add  a  viscous  damping  model  for  both  the  global  vibration  of  the  beam  and  the  local 
vibrations  of  the  delaminated  portions  of  the  beam.  Expressions  for  each  of  the  constant  coefficient  matrices  are  provided 
in  Appendix  A. 

In  the  following  analysis  it  is  assumed  that  there  is  no  uncertainty  in  the  beam  length  or  material  properties.  We 
therefore  set  L=0.24m,  h=2.25mm  E/=  75,889,600,000 GPa,  and  p  =  1234.0 kg/m.  The  damage  parameters,  Xa,Xb,a  are 
assumed  to  be  unknown  as  are  the  viscous  damping  coefficients  Ci,C2,3  (see  Eq.  (A6)  and  the  initial  global  beam  deflection 
yi(0)).  Again  we  assume  that  we  will  have  access  to  a  single,  noise  corrupted  signal  and  that  the  noise  is  additive,  iid 
Gaussian  distributed.  Thus  the  likelihood  function  is  the  same  as  was  used  in  the  previous  example  assuming  only  one 
signal  (in  this  case  the  global  motion)  can  be  recorded.  Previous  work  by  the  authors  demonstrated  that  global,  vibration- 
based  detection  was  not  possible  unless  the  measurement  was  recorded  from  near  the  delamination  site  [27].  This  stems 
from  the  fact  that  the  only  coupling  between  the  laminate  vibrations  and  the  rest  of  the  beam  is  inertial.  Thus,  for  relatively 
small  delaminations  (what  we  are  interested  in)  there  is  little  in  the  way  of  influence  with  respect  to  the  global  motion 
away  from  the  delamination  site.  The  small  influence  of  the  delamination  on  the  global  response,  even  at  the  delamination 
site,  makes  this  a  particularly  challenging  system  identification  problem.  In  what  follows  we  assume  that  we  are  able  to 
record  the  beam  response  from  the  delamination  site.  Again  we  assume  a  transient  response  to  an  initial  displacement  and 
use  N=512  sampled  points  at  a  sampling  interval  of  At  =  0.001  s.  The  variance  of  the  additive  Gaussian  noise  was  fixed 
such  that  the  signal-to-noise  ratio  was  SNi^= +20 dB. 


Fig.  4.  Actual  and  observed  impulse  response  data  (a)  followed  by  posterior  distributions  of,  respectively  (b-h),  delamination  start  point  Xa,  delamination 
end  point  Xt,  delamination  depth  a,  global  coefficient  of  viscous  damping  Ci,  local  damping  coefficient  €2,3,  initial  beam  tip  deflection  yi(0),  and  noise 
variance  cj. 
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Fig.  4  shows  the  observed  (noise  corrupted)  data  along  with  the  true  underlying  free-decay  signal.  Also  shown  are  the 
estimated  structure  parameter  distributions,  initial  condition  distributions,  and  the  distribution  of  the  noise  variance  along 
with  the  true  values  for  each  parameter.  Despite  the  very  small  influence  of  the  nonlinearity  on  the  response  the  algorithm 
is  still  able  to  identify  the  delamination  start  and  end  points  as  well  as  the  depth.  The  delamination  depth  appears  easier  to 
identify  as  evidenced  by  the  very  narrow  confidence  interval  around  the  correct  value  of  a =0.1.  The  delamination  start  and 
end  points  have  very  little  influence  on  the  global  response  as  was  illustrated  in  [27],  hence  the  confldence  intervals  for 
these  two  parameters  are  large.  Thus,  the  method  gives  us  a  clear  indication  of  the  degree  to  which  we  can  trust  our 
estimates.  In  this  case  we  can  be  fairly  certain  of  the  delamination  depth,  but  much  less  certain  about  our  ability  to 
estimate  both  the  beginning  and  end  points  of  the  delamination.  If  our  application  requires  a  tighter  confldence  interval  we 
might  need  to  try  a  larger  excitation  (allowing  the  nonlinearity  to  more  strongly  influence  the  response)  or  maybe  even 
change  to  a  local  damage  detection  method.  Regardless,  the  information  about  our  faith  in  the  ability  to  estimate  these 
parameters  is  clearly  valuable  information.  Both  the  linear  damping  Ci,  the  initial  beam  deflection  yi(0)  and  the  noise 
variance  are  easily  estimated  with  a  high  degree  of  confldence.  These  parameters  clearly  have  a  large  influence  on  the 
global  response  thus  one  might  expect  them  to  be  easily  obtained. 

It  turns  out  that  shallow  delaminations,  such  as  the  a =0.1  case  just  presented,  are  more  difficult  to  identify  than  thicker 
ones.  Consider  the  case  of  a  delamination  of  depth  a =0.2  with  the  start  and  end  points  of  the  delamination  flxed  at 
Xa=0.05m  and  Xb=0.2m  respectively.  As  in  the  previous  example  the  global  damage  parameters,  Ci,yi(0)  are  trivial  to 
identify.  The  difficult  to  identify  parameters,  Xa,Xb,a,C2,3  are  shown  in  Fig.  5.  The  variance  associated  with  the  estimates  of 
each  of  the  parameters  is  reduced  from  the  previous  case.  This  example  points  out  that  certain  combinations  of  parameters 
are  more  easily  identifled  than  others.  It  depends  to  large  extent  on  the  degree  to  which  the  estimated  model  parameters 
influence  the  observed  data.  Parameters  with  little  effect  on  the  observed  response  will  be  hard  to  identify  whereas  the 
converse  is  also  true. 

As  a  flnal  example  we  demonstrate  how  the  approach  can  be  used  to  track  damage  in  a  structure.  For  this  example, 
the  delamination  starting  point  and  depth  were  flxed  to  the  values  Xa= 0.05  m  and  a =0.2  respectively.  The  delamination 
end  point  was  slowly  varied  from  Xb=0.1  m  up  to  Xfa=0.2m.  Again,  using  only  the  noisy  free-decay  response  the  goal 
was  to  estimate  and  track  the  delamination  end-point.  Fig.  6  shows  the  progression  of  estimated  delamination 
length  along  with  the  associated  95  percent  confldence  interval.  Also  shown  are  the  “true”  values  for  Xt  used  in  generating 
the  time-series.  For  each  damage  case  we  used  a  prior  generated  from  the  previous  case.  Speciflcally,  we  selected  the 


(a)  (b) 


(c)  (d) 


Fig.  5.  Calculated  posterior  distributions  of,  respectively,  (a-d)  delamination  start  point  Xa,  delamination  end  point  x^,  delamination  depth  a,  and  local 
damping  coefficient  €2,3. 
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Delamination  end  point 

Fig.  6.  Estimated  delamination  end  point  as  a  function  of  the  actual  delamination  end  point.  The  intervals  of  confidence  were  created  as  the  central 
95  percent  of  the  sampled  values  from  p{Xb)  obtained  using  the  MCMC  algorithm. 


prior  from  a  uniform  distribution  centered  at  the  previous  mean  value  with  a  width  of  roughly  two  standard  deviations. 
For  the  first  damage  case  {Xb=0A  m)  we  also  assumed  a  uniform  prior  on  the  range  ~  U(Xa  =  0.05,0.1)  (i.e.  assume 
an  initial  delamination  length  between  0  and  0.05  m).  Again,  for  this  particular  example  the  damage  parameters  Xa,Xb,a 
have  very  little  influence  on  the  global  vibrational  response,  thus  the  confidence  intervals  tend  to  be  relatively  large. 
However  the  algorithm  is  clearly  tracking  the  progression  of  the  delamination  with  reasonable  accuracy.  By  the  time  the 
delamination  is  large  with  respect  to  the  beam  length,  the  confidence  intervals  narrow  considerably.  The  large 
nonlinearities  produce  a  more  readily  identifiable  signature  in  the  global  response,  hence  the  associated  nonlinearity 
parameters  are  more  easily  identified.  The  information  provided  in  Fig.  6  is  precisely  the  information  that  the  owner  of  a 
structure  would  be  interested  in:  the  estimated  damage  extent  and  associated  confidence  in  that  estimate.  Given  this 
information  and  the  cost  associated  with  the  structural  damage,  one  can  make  optimal  decisions  regarding  how  best  to 
maintain  that  structure. 


6.  Conclusions 

This  work  has  presented  an  approach  for  identifying  nonlinear,  multi-degree-of-freedom  structures  using  only  observed 
free-decay  response  data.  The  approach  appears  to  work  well  for  limited,  noise  corrupted  observations  that  are  easily 
obtainable  in  experiment.  This  makes  the  approach  attractive  from  a  practical  standpoint  with  the  main  drawback  being 
the  computational  effort  required  to  build  the  stationary  Markov  chains.  At  each  step  in  building  process  the  practitioner  is 
required  to  numerically  integrate  the  model  forward  in  time  for  each  parameter  being  identified.  Nonetheless,  the 
approach  is  effective  at  identifying  both  linear  and  nonlinear  model  parameter  distributions  for  structural  models.  This  is 
precisely  the  information  that  is  of  interest  in  damage  detection  problems.  Here  we  have  show  the  approach  can  be  used  to 
estimate  and  track  the  length  of  delamination  in  a  composite  beam  using  only  noisy,  free  decay  response  data.  When 
combined  with  knowledge  of  the  costs  associated  with  structural  damage,  this  approach  provides  the  information 
necessary  for  making  optimal  maintenance  decisions. 
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Appendix  A 


Eq.  (10)  was  derived  using  energy  methods.  The  constant  coefficient  matrices  associated  with  these  governing 
equations  are  given  here.  The  linear  stiffness  matrix  [Kl]  is  given  by 

/  hh^(-4L3+3(a-l)a(Xb-Xa)^)E 


[Kl]  = 


V 


1616 

0 

0 


0 

a^bh^n^E 

6(Xb-Xaf 

0 


0 


(Al) 


The  two  matrices  involving  quadratic  terms  are  defined  as 
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V 
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3(1  -a)abh^  n^(2L-Xa  -Xt,)E 
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f  3(\-a)abh^n^(2L-Xa-Xti)E  3(1-a)abh^n^(2L-Xa-Xi,)E\ 
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0 

Vo 
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0 
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(A2) 


The  matrix  multiplying  the  cubic  term  is  given  by 


[/^c]  = 


/O  0 

Q  abhn^E 


SiXb-Xar 


0 


0 


0 

0 

(\-a)bh7i'^E 
8(X6-Xo)^  J 


(A3) 


V 

The  mass  matrix  is  slightly  more  complicated.  Defining 

j^^^bh(Xi,-x„)[(x„+Xi,)p((-3  +  7i^)xl+6XaXi,+(-3  +  n^)xl)-2L((-3  +  2n^)xl  +  2(3  +  7i^)x„Xi,+(-3  +  2n^)xl)] 


the  mass  matrix  is  given  by 


[M]  = 


140 


bhLp 


aoc 


(l-a)a 


-abhpiXb-Xa) 


0 


(l-a)a 

0 


-('l-a)bhp(Xt,-Xa)  j 


(A4) 


(A5) 


The  damping  for  both  the  global  beam  motion  and  laminates  is  assumed  to  follow  a  viscous  model  such  that  we  have 

/Cl  0  0  \ 


[C]=  0  C23  0 

I  0  0  C23 


(A6) 
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